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This is the Final Technical Report for the NASA Cooperative 
Agreement NCC 5-14 "Application of Laser Ranging and VLB! Data 
to a Study of Plate Tectonic Driving Forces" for the funding 
period 1 May 1979 through 30 April 1980. The broad goals of our 
research under this Cooperative Agreement have been to investigate 
the conditions under which changes in plate driving or resistive 
forces associated with plate boundary earthquakes will be 
measurable with laser ranging or VLBI and to pinpoint those 
aspects of plate forces that can be characterized by such 
measurements. 

This Report is divided into three parts, each describing 
a specific research task completed during the past year and 
contributing to the broad goals of the project: (1) analytic 

solutions for two-dimensional stress diffusion in a viscoelastic 
plate following earthquake faulting on a finite fault; 

(2) finite-element solutions for three-dimensional stress 
diffusion in a viscoelastic Earth following earthquake faulting; 
and (3) quantitative constraints from modeling of global 
intraplate stress on the magnitude of deviatoric stress in the 
lithosphere 


PART 1; ANALYTICAL SOLUTION OF TWO DIMENSIONAL TIME 
DEPENDENT STRESS PROPAGATION 

1. Intrcducticn T " 


Most of the eiirthquakes today occur along plate boun- 

' ■■ '■ I 

daries. Strain energy accumulates and releases at these 
bou^ndarios /deforms surfaces/ and - causes dammage. The 
understanding of earthquake mechanisms is bs^qcraing more 
iinportant in order to forecast and, further/ to predict 
earthquakes. . .-V:. ■ . ' 


The rcigrati.cn of laege eartligliakes was first studied 
fay I djogi ( 156 8a / I 968b) . Delsemme and Smith (1979) -have 
qua|ntitatively analy.z.ed the migraticn pattern: along many 

I ^ ^ I 

plate Loundaries aroun^d the world; . periodic trends in 

■ : ,i ■ i I 

space and time ate obtained. A qualitative mechanism con- 
cerning a possibly deep-seated migration front of stress 

_ i' ' ' ' 

triggering earthquakes along its propagating direction 

(Scholz/ 1977) is jthe general explanation to the migration 

r ' cm.!' . 

phenomenal However/ the relaticiiship between large earth-:: 
quakes at greater tim e-distance intervals is less clear. 
Recent ly/" the nigration of crustal deformation has been 

i 

discovered in many tectonic areas ( Kasahara / 1 979) . 
Ijispersion and dissipafcion of the deforiration waveform are 
also riotod as chacacl erist ics . ^xtrapolatinj the ftijra- 
tion FMth Lack toward the ocean , the de format ion frdnt 
soems to have o rig ir.at^i from the vicinity of thji trerich. 
It is supposed that cither a repeated irrog jla r a scisn ic 
plate Ecticn generates the def etna tion e VentS/ or that it 
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ce.sults fcon a periodic seismic slip at a plate boundary. 


I 

I'- r 


Qudiiti tatfve study has teen carEied out by Yang (1979). 

Analytical sc luticn a.boct the^ stress (or I displacement) 

tesporse ini! time’ a bd: space due to earthquake disturbances 
at the plate boundary i s still hedessary to obtain physi- 
cal insights of the! inigrat ion Lpf ccess.__^^ 

I n^^^ p^ reports, vie have studied the effect of 

earthquakes on plate mcticn using analytical one- 
dimersional viscoelastic and tw o-dimensicnal elastic 

' ■ i i ■ 

models. The itcdols inccrpccated the interactions between 
the elastic; lithosphere and the as thenosphere. Rijth 
poriodid earthquakes the one dimensicnal study show signi- 
ficant fluctuations cf displacement field at large dis- 
tances from the plate boundary. The two dimensional model 
consists of two finite elastic rectangular plates ,cne 
overlying the ether," with different thicknesses as shown 
in figure 1. Applying a unit diSplaceraent at the plate 
boundary suggests that the instantaneous elastic dexorma- 
tions are si jn if icant wi thin a distance equivalent to the 
Icn'gth 'Of ■;'the''''fa"ulb",.;ruptuce^ . 

In this report, we extend our two dimensional Glastic 
stuiiy to a visccelastbc" moi;eJL using the corcedpoadonce 
principle (Christiansen , 1 97 1 ) . A jlaxwollian viscoelastic 
a;;thi:r.osphere is assumed and t ime- dependent ;io forma ti c ns 


in the lithcsnhere are calculated. 
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2. Tvo Dimensional linie-Depe nde nt Solution of Stress Pro- 
pagation , , 

i'_ ^ ‘ 

^ c-:- '■ r.J I 

Tine aependent scluticns of - ; horiaontal stress and 


strain propagation i,n lithosphere candle obtained using 

correspondence principle (Chf istians'en, 197 1)* Figure 1 

• i ' ' 

shows a finite two diraensicnal plate wifh length L, width 
2C, thickness Hi overlyincj a second layer of the same size 
and thickness H2. We choOjSe a simple' iiax well ian Viscoe- 


lastic behavior in shear for the 


asthengsph^re. 


and ah 


elastic response in the lithesphere. The hatched area 
ropre'sents the finite rupture zone at the plate boundary 

with total width 2D. The x-coordinate is normal to the 

• i ' ■ ' 

plate bcundary and the y-cccrdinate""]is pdrall;el to the 
plate boundary. The originj is Iccated at the center of 


the "fault". V.’e are calculating the stress and displace- 
inent at specific distances frem the origin in both x and y 


di rectiori s. 


Starting with the eguatiens of notion in the two 
dinicnsicnal elastic case ; 

2 MT-w ) luv) a^v/ {^x . 3y) 

-2/( ( 1 -v'’ ) u/ (E • FI • H2) (ba) 

2 d'v/jy"t(1-K') 9^'/ u/(jx-^y) (1b) 

-2/<( l-vV) V/ (E -HI •Ii2) 

vi.uro u, V arc disp) la cement field in X, y ilifection 
L-Oor'.-ct ivuly* V is the Pcissor.'s catiov S is the Young's 


modulus of the lithosphere, and yu is the shear modulus of 
the as thono sphere. 

Tahe Laplace transform of- (la) and (1b) and ncn- 
dinensicnalize - 

- I ! ■ ■ 

u, V by the in itial displacemen t at rupture zone, 

X, y by li1, the thickness of the lithosphere, and 
t by T, the relaxa tion time of the asthenosphere, 

wc obtain: — ■ , — ' 


2J'* (IpV) (1+V) 3^V/ (3 x-3Y) 

= 2 (Hlyup/ (Si- H2)-:, (l-v* ) T- s -U/(T-ai"1) 

2 9 v / 9 y ^+ (1-v') 9*v/Jx*t+r^(i+w) ^\/(9 x-3y) 
= 2 (Hl-yx) / (B‘H2) (1-V*') Ivs^^;V/(T-s+1) 


(2 a) 


(2 b) 


viiere X and Y are tlie non-dimensional xt and y j d and V are; 
tl'.e non-di mensicnal u and v in Laplace domain, T is the 
relaxation time of the asthenosphere, and s is the Laplace 
r-a rare tor. 


i \ T 


The above two aquati.ons ate similar to those in our' 
two diiiensional elastic case except that there is a time 
deponcerit factor Ts/(Ts+1) at the right hand side of the 
e» 7 u.i cions. 


;;e can fuLtlioc simplify the e-juntions (2a) and (2b) 
inter . 


I 


i 




! 

i 

I 

I 
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=k‘^/2 

(3) 


(4) 


vhere 7*" is the laplacian operatcr in two dinonsicn and 

f- = 3 u/ax+ 9 v/ay 

j ■ ■ ' 

= 9 U/, 3 Y-aV/ 3 X - 

r::i" 

k"" (s)=2(/^'H1)/(E-H2) n-»»^)Ts/{Ts + 1) 

A set of solutions that satisfies equations (3) and 
(4) is 

Wfv.J) = £ 

+ Z-SVns<;i 
»M =0 y 


)J 

0. cv< -ivX U / a ,>^^4. L) > 

^(Aje -^e ) - -^r^jr^e.e )j 


Oo t rj aV -<0^ b bx -bvi 


iv. , / -cO( U ox . ) 

+ s; *>'Ae >[ 

tM=o (/I CC~<JL ” ■' 


bx ~bX 
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The coefficients A1,ft2, are functions of m. 

Notice that The displacement fields decay exifonentially 
with X'and' the decay constants, a and b, are tiiac- 
dcpc-ndert. ■ 

Boundary conditicns / - 

1' ---• 

Great earthquakes and their aftershocks onT-y occur 
alonij-^ a saiall portion of the toundary. To simulate a fin- 
ite rupture zone, we apply a unit displacement in the 
hatched area from y=-D to y=D throughout the” thickness Hi 
(See figure 1) . remainfdcr of this boundary is fixed 

in the x direct icBi ; Figure 2 gites 'a top view of the 

model arid a brief summary of the boundary 
Motion in the y direction j,s iallowed at the trehTh, The 
rillge side is asisured tc be stress free. Along the other 
two bo uuvdaries,' we '"confine the motion to the y direction 
and ignore friction. 

- ' -- 4 '' ■■ 

; ! ■ ' - 

Ha themat ica 1 ly, these conditions are: - ' ' 


U (C, Y) =D/C+ 2 sin go? ( oC-Y) / { ck C) 

d' X (L,Y)= ^Xy (L,Y) =0 



Using' these boundary conditions and eguaticns (3) y 
(4), we obtain tise-dependen t coefficients with the s-kme 
forn as those in two diraensionai'ela case . Inverting 
to ! tl.c tir;e doir. ain from the Laplace dcniain with the Fast 
Fourier transform, simple numerical solution results 
(Dubner and Abate, 1S68) . 


Pe suits 


In most of. the figures shewn below, ve use two kinds 
of units, dimensionless unit and dimensional unit. For 
dimensional unit ve refer to 


L=5000 k m 
C = 4000k, ffi 
HI = 100km 
I!2 = e0km 
D=300kra 

/<= . 7 e 1 1 d y ne /c m 
F= 1 . 7 e 1 2 dyne/ cm 
1 = 1 el sec 
30 1 9 poise 


!'.’e assume 1. 6 met ers horiz c ntal ' disp lacemff nt normal 
to the rupture zone to simulate earthguake slip. 


The u displacoment parallel to the X-axis at dis- 
tances of 4 00 , 00 0, an d 1 2 00 k ra frem the trench are shewn in 

. / . ■ ^ ' . . ; ' ■■■' - , . ■ ^ , ' ' , , ,, 

fijure 3. At time 0, the .ilstur Lance at plate bcundary 
effects only adjacent regions of the plate. As tine 
iiic r(>!ast.-s, it graiiually propagates dcreper into the plate. 
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The velocity as a function c£ time alony the x-axis 
is shewn in figure 4. The effect of the asthenospheric 
str[es& relaxaticn is obviohs. The shape of the velocity 

curve is similar to those in cur one dimensicnal moiel,: 

I ■ ■ ' I . i . ^ 

w hich assumes an infinite taijlt; however, the magnitude of 

the maxirau)! ; velocity is smaller than the one dimensional 

’ . - I I 

case. for example, the maximum velocity at 400 km from 

the fault is about 251 small. This suggests that fault 
length is one impertaat factor for propagation velocity. ~ 

. ! I : 

riguco 5 shows the v displacement along the positive 
y axis at | time 0,5, 10,20, 4 0. The maximum v-displaceaeiit - 

1 ^ I 1 j I 

always occurs at the end of the rupture zone. The dis- 

placemontj field decays e xcc^entially -,wlth time. The 

relaxation effect is not obvious in this direction because 
our initial displacement is essentially perpendicular to 
this direction. 

At the tinie of the earthquake very high“ shear ^ 
stresses (about 110 bars) occur at the end of the rupture 
zone as shown in figure 6. The negative sign represents 
the orientaticn of the shear. These are dependent on the 
Lo.ir.ni'.ry cop.il itiens . 

The time depen-Ient solution of stress TT ' at . di.nen- 
siot.l(-ss distances 0 , 4, 0 along t he x ax is is slicwn in fig-i h 
ure 7. “or a 1-6 meter horizontal slip at the botnd-.iry. 




the istress at x=0 is about 22.5 bars, which is consistent 
with finite; element calculations (Yang, 1979). The stess 
decays exponentially with distance for s^ort time inter- 
vals after the earthguake disturbance. At several hun- 
dreds of kilometer away frcir. the boundary , the relaxation 
effect is again significant. 

Figure 3 iis the distrituticn of stress CT x along the 
y-axis. For short times, the stress O^x changes sign 
across the Iboundary of the rupture xcne^ at y=300km. A 
sudden slip in the rupture zone inSeases^ t st res s level 

outside of the zone parallel to the plate boundary. -his 
slight increase in stress may trigger earthquakes in the 
neighboring area along the plate boundary if the tectonic 
stress level is already large. 

3. Summary 

Our two dimensional analytical study of stress propa- 
gation has been completed. Although this is a simple two 

i ■ ' 

layer mcdel, it gives results that retain the significant 
features found in more complicated tin itje element ircdels. 
It else provides seme pliysicai insight such as the time- 
aependent oxpcncntial dec^l y constant of the displucement 
a:ui stress fields and allow us to understand the process 
ill some detail. liov.ever, it is not easily generalized to 
GpoGific rogtens. Cur future study will conceritrate on 



the •devolopmenti of an ‘infinite element' based 
analytical sclutions to simplify three dimensional 


cn cue 
finite 


element modelling. 
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FIGURE 1 

Finite twc-di nensicnal plate; inodel with “ length^^^^L 

I i , . 

vioth 2C, ' thickness _of the lithospJiere Hi, and the 

cisthen csphreic thickness H2. The hatched . a tea 

1 r ' * " ' „■ 

repcesents the finite repture 2 one at the plate boun- 
dary with width 2D. The x-coordinate is iicrnial to 
the plate boundary and the y-coordinate is parallel 
to the plate boundaryi. The origin is located at- the 
center of the "fault”. Se apply a unit displacement 
throughout the thickness of the lithcsphete in the 
rupture v:cne tc simulate an earthquake. 




' - figure^ 2 - — : 

This is the top view of the moflel. A unit- di^place- 
Ront in , the rupture zone norinal to the trench is 
applied to simulate earthguahe slip. Ihe rest of 
this fcoundary is ^fixed to the x-direction; 
motion in the y-c^ection is allowed. ^Th^ i^idige side 
is assumed to be Stfesjs free. Along the other two 
boundaries, we confine the motion to the y direction 

and ignore friction. 


JJrtX 



BOUNDARY CONDITIONS : ' 

10 UNitl X-DISPLACEMENT IN THE RUPTURE ZONE 
2D RIDGE IS STR;ESS FREE , ; ' • '^ ! | 

3D ZERO Y-DISPLACENENT ON THE OTHER“TWO“BOyNDARIES- n 

: , ii 

" 

Figure 2 . — ^ 
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The u- displace!no, 7 t parallel to the x axis. Distaapes 
are rion!“d|Lnienai 9 :n:Tali 2 €d fcy jt hie t-hicknegs | of , the 
llthoaphere . The number by the curve is the rdiraen- 
siorilessl -time. At time 0, the disturbance at plate 

Lounjdary effects only adjacent regions. As time 

increases, it gradually propagates deeper into the 
plate. . ; 


! FIGURE 4 ' 

The velocity (time derivative of the u-displacenent) 
as a function of - time- alo n y the x-ax„is at jc=4,8,Jl2| 
T h'e e fif ec 1 1 of the asthenos p hef iff stress relaxation is 
olivious. The ma xim um velocity at x=4 is about 25S 
smaller than the corresponding one dimensional case. 







FIGURE 5 


riie v-displaceisent along the pos±t ive y dlrocticD at 

1 ■ ■ — • 

0,5,10^20^40. The aaxifiuiiTi displacement occurs 

'T* ' • - ' 

at the end of the rupture zsRe. Displacedent ; decays 
rapidly at short time., The relaxation effect is not 
obvious in this direction because of the direction of 
initial conditions. 


■ , :;r^' FIGURE •6 

Snear stress along the positive y 1 axis. The dinen- 
siciKil units is obtained using the parameters speci- 
fied in the text.. Very high shear stress occur-- at 
the end of the rupture zone. The negative .dign 
reprc?sonts the crientaiton cf the shear. 
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FIGURE. 7 

liie stress d^x at diinensionless distances 0,4/8 alcng 
the axis. For a 1.6 meter horizontal slip at the 
uaxiirura stress ■ right behind the rupture 
sor.e, is about 22.5 tars. At greater distances/ the 
stress relaxation is significant. 


FIGURE 8 

ristributicn- of stress dx along the y axis. For 
snort timos/ the stress changes sign across, the boun- 
dary of the rupture zone. A slip in the negative x 
direction increases the ccmpressionil stress outside 
of the rupture zone parallel to the i plate boundary. 
This ffight trigger €^rt]^guak6s where the tectonic 
stress is alrGady high.. 
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PART 2; TIME DEPENDENT DEFORMATION AND STRESS RELAXATION 

AFTER STRIKE SLIP EARTHQUAKES 

1. INTRODUCTION i 

Great damage has been caused by shallow strike slip 

" ■ 1 

earthquakes along plate boundaries in various parts of the 
world. The mechanism of these earthquakes has long interested 

i 

\ I' 

seismologxsts. The study of geodetic measurements of the 

I i , : — ■ !■ 

1906 San Francisco earthquake led to the formulation of 
'll .i ' ! I : 

elastic rebound theory (Reid, 1910) r much of which has 

; i ’ ■ ^ 

remained a basic tenet on earthquake mechanism. The 

continuous study of accumulation, release and relaxation 

of stresses near the fault zone has provided a more detailed 

mechanism of strike slip earthquakes (e.g.: Chinnery; 1961; 

Scholz and Fitch, 1969; Turcotte and Spence, 1974; Savage,; 

1975; Thatcher, 19 75a, b; Budiansky and Amazigo, 1976; Rundle 

and Jackson, 1977a, b; Savage andj Prescott, ,1978; Savage ,| 1979; 

Thatcher, 1979; Turcpotte et al., 1979). Much of the study 

used geodetic measurements near the fault zone. In particular, 

' i 

static elasticity and dislocation theory have often been applied 

to the study. -i™ . i \ ;;; r 

Stress, accumulation, release and relaxation are time. 

I I : 

dependent phenomena. This is evident from geodetic data, the: 
migration behavior of , earthquakes and asthenosphcric viscosity. 
In recent years, there have been intensive geodetic and creep 

measurements in the San Andreas fault Zone, and ultra precise 

_ 1 ! . 

space technology has been applied to geodetic measurements 
(e.g. Niell et al., 1979; Smith et al. , 1979) .■ -Accurate 


2 . 





data will be available in the near future on the time 
dependence of crustal deformation. A detailed, three- 
dimensional time dependent model will be necessary for 
interpretation of such data. On the other hand, earthquake; 
migiration phenomenon have been observed along plate bounda,ries, 

I ^ ' 

most noticeably along the North Anatolian fault (e.g. Mogi 

i' „ 

1968; Allen, 1969; Dewey, 1976; Toksdz et al., 1979). 

Explaining this time dependent phenomenon also calls for 
time dependent models. 

■We calculated the long term time dependent response of 

, r::: 

a s^t of models of strike slip events. The effect of 
relaxaltion is isolated for the calculation. Most of the 

I 

previous attempts to model time dependent tectonic phenomenon 
after earthquakes used two-dimensional models (e.g., Nur 
and MAVkOF 1974; Bischke, 1974; Smith, 1974; Savage and 
Prescott, 1978; Thatcher and Rundle, 1979; Melosh and 

'll '1 : 

Raefsky, 1979) or simple layer and half space solutions 

i : . . • ' “ 

(Roseninan and Singh, 1973a, b; Barker, 1976; Rundle and 
Jackson,. 197 7a, b; Cohen, 1979; Cohen and Cook, 1979; 

Lehner et al., 19 79). Two-dimensional models_ass;ume 
an infinite long fault, and the effects in the region beyond 
the fault tip cannot be described. However, in this pape r- -i v 
we show that there are significant effects in the region 
beyond the -faifit tip- for strike slip events. Laterally 
homogenebus models that assume no lateral heterogeneities ^ 

across the fault zone oversimplify the near-source probiem. 
Most data indicate the presence of lateral heterogeneities 
^ near . the-:<faultv' ^ 
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In this paper, we present time dependent calculations 

; i ' ■ 

for finite strike-slip faults in laterally heterogeneous 
media. We use the three-dimensional fiinite element method 

to model sttike slip events. i The forward problem is set up 

i - . ‘ ' I 

to predict time dependeht deformation and stress after a strike 

slip event for years and tens of years. The models are 

( 

chosen for representative earthquakes to show the characteristic 
time dependent features of strike slip events'. The boundaries 
of inhoniogeneities in the models are kept geometrically simple. 
The model results indicate that geodetic measurements after 
an event may provide information on rheological properties 
near the , fault zone which are vitally related to earthquake 
occurrence. 


II. MODEL DESCRIPTIONS 

I : ' ' ' ' 

Fault Models 

We present the computation for two classes of models, one 
for a great earthquake, and the other for a moderate size 
earthquake. Due to uncertainty in the viscosity structure , 
several sets of viscosity values for the same fault model 
are used to show a range of possible results. We are 
interested in the general behavior of f elakation. No attempt 
is made ' to model a specific region in detail. However, the 
fault displacement and dimension for great earthquake models 
(Gl, G2 and G3) are comparable to those of the 1906 San 
Francisco earthquake or the 1939 North Anatolian earthquake. 
The dimension for fault models Ml and M2 are appropriate 


4 . 



for a magnitude 5* 5-6.0 earthquake , such as the Coyote Lake 
earthquake of 1979 (M = 5.7). 

We model the strike slip earthquake as a sudden slip 

on I two sides of a rectangular, vertical fault surface. The 

offset vaiLues are prescribed, and the subsequent displacements 

and relaxation are computed. We assume th^^^^ the slip on 

fault stops after the step function slip. This probably 

will happen, if deviatorib tectonic stresses are relieved 

near the fault zone jand friction again takes over, and the 

recfion deforms in a coherent manner. However, there may be 

i |: _ ^ - 

cases where this condition is violated. Slow afterslip is 

sometimes observed after an earthquake (e. g. , Burford, 1972; 

Buckham et al., 19777 Goppersmithj e!t al., 1979). If 

afterslip : does not last long relative to relaxation time 

(years) tllie I long term ef fect will be similar to a step 

function. 

The mode l^region for models Gl, 02- and G 3 is 2740 km 
long, 2320 Jen wide and 700 km deep.; The fault is a rectangle 
350 km long and extends vertically from 0 to 40 km depth. 

The relative fault offset is 5 meters sl[rike slip, i except 
that nearir the edges of the fault area^ it tapers off . ‘ The 
offset tapers off linearly from 5 meters at 20 km depth to 
zero at 40 km depthr-'and it tap linearly to zero at 

35 km from the tips along the strxke direction as shown in 
Figure 1. 

The model region for models Ml and M2 is 196 km long, 

169 km wide and 80 km deep. The fault is a vertical 
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rectangle 20 km long and 12 km deep. The model earthquake 

I,.'"' ' , . 

has no surf^e rupture. It has relative offset of 30 cm 
from 3 to 9 km depth, and the offset tapers off linearly 

to zero from 9 to 12 km depth and 3 4o 0 km depth. It also 

I 

tapers off linearly to zero at 2.5 km from tips along the 
strike direction (Figure 2) . 

Material Properties 

There is controversy over the laws governing creep 
behavior of earth material (VIeertman, 1978) . Linear Newtonian 
behavior of the mantle fits the post-glacial rebound data, 
but laboratory rock mechanics experiments show non-linear 
creep behavior. For calculation of perturbation caused by 
the earthquake, we choose the linear viscoelastic models, 
simply because viscoelastic material contains the essential 
elements of time dependent relaxation phenomena, while the 
assumption of linearity greatly simplifies the physical 
picture, as we can separate the effect of perturbations 
caused by the earthquake. 

The material in the model is assumed to be elastic in 
bulk and maxv/ellian viscoelastic in distortion. The short term 
elastic constants of earth vary slowly in space (Hadden and 
Bullen, 1969) . On the other hand, the viscosity value changes 
by orders of magnitude from lithosphere to asthenosphere 
(Cathles, 1975; Peltier and Andrews, 1976) . The contribution 


to viscoelastic relaxation due to changes in elastic parameters 
is therefore relatively unimportant. In order that we do not 
unnece„ssaiLily complicate the physical picture, we assume in - 
all the models that the elastic parameters are uniform in u 
the region, with bulk modulus 1.3 x 10^^ dyne/cm^ and Poisson's 
ratio 0.25. ! The different model results will be due to 
different viscosity structures in the models. 

The viscosity of the earth is not a well constrained 


quantity, especially near a tectonically active zone such 


as a transform fault. The lithosphere in general can 


wx 


deviat 


oric stresses for aj million years or longer. 

! 


The thickness of thei elastic litHos|phere i-n‘ continental 


regions is probably greater fhan 50 k^. The as theino sphere 
has' a .visOoiity vailue Orders 'of -rmagnitude smaller Ithan the 


lithosphere. From analyses of post-glacial rebound data, 

the low viscosity layer beneath the lithosphere can be on 

20 ■- 

the order of 10 poise (Cathles,_1975 ; Peltier and Andrews, 


1976). The 16w viscosity implies that deviatoric stresses 
cannot be sustained there for a time scale larger than several 


years to decades. Ur-j ■ ^'77' ::^r ■ 

Near I an active transform fault, where only shallow 
earthquakes are observed suggests lateral heterogeneities. 
For examples, in the San Andreas fault zone earthquakes are 
usually shallower than 10 to 15 km, indicating that stress 
is being relieved anelastically below this depth. This 
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thickness of the "seismo-genic” layer 'is simply ioo small 

compared to the generally accepted lithospheric thickness, 

I 

suggesting that a low viscosity zone may exist under the , 
fault.! Lachenbruch and Sass (1^73, 1979) found that thie 

1— z. _ ' ' ; ! ! ■ i 

San Andreas fault system is contained in a broad zone df 

i ! ! ; ; 

high heat flow anorftaly. They concl4ded that the thin 
seismo-genic layer ! is more" brittle than the layer beneath'^^^ 

!i ■ , ■ ' i : . 

it, implying the possibility of shallow low! viscosity zone 
there. Theoretically, the shearing 'motion should also cause 
temperature and viscosity-structure to vary laterally (Yuan 
et al*, 1978). Three dimensional inhomogeneities in elastic 
parameters are also observed in tectonically active regions 
(Zandt, 1978) although it is difficult to estimate viscosity 
structures from elastic parameters. 

There is little information on the viscosity in a 
fault zone. Budiansky and Amazigo“ (1978) estimate the 

M 21 

"effective viscosity" of lithosphere to be 10 poise in 
California. Nur and Mavko (1974) and Smith (1974) analyzed 

! i 

the vertical deformation of the 1946 Nankaido earthquake and 
conclude the viscosity value below the elastic lithosphere Is 

on the order of 10^^ to 10^^ poise. Thatcher and Rundle (1979) 

20 ■ ■ ■ 

found that a viscosity value of about 5 x 10 poise in the 
asthenosphere near Japan fits geodetic measurements. It 
is plausible that the shearing motion and higher temperature 
in the fault zone make the viscosity lower than in its 
adjacent regions. We will carry out calculations for a 
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range of viscosity values. 

Great earthquake model G1 is the control model; a j 

layered structure is assumed- The lithosphere extends to [ 

25 ' 

80 km depth; it is given a viscosity value of 10 poise i 

from 0 to 40 km, and 10 poise from 40 to 80 km depth. | 

A low viscosity layer extends from 80 to 180 km depth, 

with a viscosity value of 10^® poise. Below 180 km, the 

22 

mantle viscosity is assumed to be 10 poise as shown in 
Figure 3 . 

' i 

More realistic models incorporate lateral heterogeneities 
across the fault. In model G2, a low viscosity zone rising ^ , 

to 20 km below the surface and extending 140 km on each side 
of the fault is assumed. The viscosity value is assumed to 

20 I 

be 10 poise, the same as in the low viscosity asthenosphere 
extending from 80 to 180 km depth. Calculations show that 
differences in resulting time dependent deformation and stress 
relaxation between model G1 and G2 are significant. Model G3 i 
has properties intermediate between model 1 and 2, with the ■ 
low viscosity layer extending to 40 km depth. For '“moderate ' 
earthquake models Ml and M2 we assume that under the fault 

the low viscosity zone extends; to shallower depths. Far 

i : ' I . 

away from the fault, the lithosphere ; is 80 km thick. Both 
models assume* that from 12 km to 46 km depth there is a 
low viscosity zone 40 km wide on each side of the fault; 
from 46 km to 80 km depth, it is 70 km wide as shown in 
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jj 


(C 


Figure 4. Fast relaxation time is assumed for model Ml 

-- ig r-.-:, -i;:.-:-:- ---■ 

(vxscosity 10 poise) to establish the maximum possible 

effect of relaxation. Model M2 assumes the low viscosity 

20 . — 
value to be 10 poise. 

Computational Scheme 

We use three-dimensional time dependent finite element 
models to calculate time dependent motions following an 
earthquake. We combine the frontal solution technique 
(Ironsi, 1970) to the unified time stepping approach , 

(Zienkiewicz and Cormeau, 1974) for a versatile and efficient 
solution scheme. The calculation scheme is described in 
Appendix A. 


Due to syii(mi^try of a: vertical strike slip earthquake, 

only a quarter of the region is needed to be modelled in 

• : . ^ ' 1 

the numerical ScHeme. All the models in this study use 720 

elements and 2982 degrees of freedom. The grid structure of 

moderate earthquake models is a scaled grid of great 

earthquake models;. Figure 5 is the top view of the grid 

structure 6^ both classes of models. The three-dimensional 

r i i ’ ' 

model is made up of^ seven identical grid surfaces. The 
element is jin 8 node, 24 degrees of . freedom hexahedron 
with; 8 ; gaussian integfation stations. Elements with 27 
integration stations have been used on some models. The 
results are nearly identical with those of 8 integration 
stations . 
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since only a quarter of the region is used in the 
■ ■ ■ ' i i. 

numerical scheme, the boundary conditions must simulate 

those of the complete strike slip fault. Thei fault area 

is contained in the boundary plane x - 0 (Figure 5). The 

boundary plane y = 0 is a symmetry surface bisecting the 

physical fault model. On plane y = 0 the displacement in 

X direction is constrained; on plane x = 0 the displacement 

in y direction is constrained, except on the fault surface 

half the fault offset value is prescribed. The effect 

caused by the event decreases with distance from the fault. 

We choose a region large enough such that the artificial 

external boundary condition does not influence the behavior 

of the region v/e are interested in. We used rigid and 

free boundary conditions on the sides of the region. We 

present results only for those elements where resulting 

stresses differ by less than a few per cent for these two 

extreme cases. We used 12 time steps for models Gl, G2 i 

and G3 to calculate time dependent values for up to 49 

years after the event. For thin lithosphere and fast 

relaxing models Ml and M2, we used 14 time steps for 

9.5 years after the event. 

III. Model Results 

Great Earthquake Models 1 7 ; i 

The results of time dependent deformation and stress 
from different models are presented and compared in this 
section. Model Gl is the control model with a laterally 
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homogeneous layered structure. The horizontal displacements 
on free surface at selected locations for model G1 are 
shown in Figure 6a. The four figures give the horizontal 
displacements at time = 0, 9.5, 25 and 49 years after event, 
respectively. The initial (time = 0) pattern is typical of a strike 
slip fault in elastic media. Afterwards the displacement 
gradually increases with time. The initial displacements 
decrease very fast with distance from the fault. The 
relaxation process [spreads the deformation outward from 

i ; 

the fault zone . i 

In contrast with model Gl, model G2 has a low viscosity 
zone extending to shallow depth near fault. The result of 
horizontal displacement on free surface at selected locations 
is given in Figure 6b. The instantaneous response is the 
same as for model Gl, since they have the same elastic 
parameters. However, the magnitude of time dependent 
displacement is much larger and concentrated near the fault 
zone (where the low viscosity zone is shallow) . The time 
dependent effects can be seen more clearly in Figures 7a 
and 7b, where the displacements albng the line perpendicular 
to the center of fault (x axis in Figure 1) are shown. Near 

the fault zone the tiine dependent Reformation is in general 

: 

small compared to the instantaneous response for model Gl, 
while the shallow low viscosity zone in model G2 significantly 
increases the magnitude of time dependent displacement. 
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The contours of vertical displacement on free surface 
immediately af±er the event are shown in Figure 8a. Also 

! 

shown (Fig, 8b) are the shear stress contours. The results 
are what we expect from a shallow strike slip event. In 
the vejrtical displacement figure there is subsidence in a 
broad area in the upper right quadrant, except near the 
fault tip. Our grid resolution cannot resolve the very small 
uplift zone near the fault tip, but the overall pattern is 
very similar to that of the analytic half space solution of 


Chinnery (1961). Subsequent time dependences show 
characteristic differences between layered and lateral 


inhomogeinepus models. The contours of vertical displacement 
change after the event (total displacement minus instantaneous 
response) at 5 and 25 years after the event for model Gl and 
are given in Figures 9a and 9b. The results of model G2 are 
given in Figures 10a and 10b.. On the upper right quadrant, 
for model Gl the time dependent vertical movement is continous 
subsidence near the fault zone, and uplift away from it. This 
result is expected if the horizontal displacement is 
spreading out from the near fault region. The magnitude of 
this vertical movement is on the order of several centimeters. 
In contrast, the results for moidel G2L is continuous uplift 
throughout the upper right quadrant. The physical •'reason 
for this difference can be seen from the horizontal 
displacement plots (Figures 6a and 6b>. In the upper right 
quadrant in the figure, the "flow pattern" of the strike 

i 

slip event is such that material 
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enters the n^ar fault zone parallel to the strike direction, 
but leaves the fault zone in the direction 45 degrees from 
the strike direction at large distances. The low viscosity 
channel near the fault zone makes it easier for material to 
enter the fault zone; the thicker lithosphere beyond the ' 
low viscosity channel forms a barrier for material to fan out, 

thus material piles up near the fault zone and a bulge 

■ ■■ ■ ’ i ’ ^ i i 

results. The uplift in model G2 is about 15 cm in 25 years. 
This effect is measurable by geodetic means, and could serve 
as a toor for investigating the viscosity structure of the 
fault zone. 

Another interesting phenomenon is the time dependent 
character of the stress relaxation for these two different 

: i ... - ^ ^ 

models. The instantaneQus perturbation of horizontal shear 
stress component at 10 km depth caused by the strike 

slip event is shown in Figure 8b, and at 25 years later 

is shown in Figures 11a and 11b for models G1 and G2 
respectively. Before the event, due to relative plate ’ 

motion the stress component a-v^r should prevail near the 

' ; - . ■■ , ““ 

fault zone. The earthquake relieves the prestress along the 
fault, but reenforces the prestress in front of the fault 
tip. ; More detailed timei dependency can be seen in Figures 
12a and 12b, inhere a^y Vs.; time at selected positions is shown 
It can be seen that along the side of the fault, a is 

negative, implying the prestress is relieved. However, 
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viscoelastic relaxation in general accelerates the stress 
recovery. In front of the fault tip, is positive, and 

. I 

the prestress is reenforced. This result of reenforcement 

^ I 

of prestress has been considered to be the cause of secondary 
faulting (Chinnery, 1966) or creep and aftershocks in front 
of fault tip (Scholz et al. , 1969). The earthquake may in 
time trigger subsequent events along the same fault. 

However, for the uniform thick lithosphere model Gl, the 
time dependent stress is small compared to the instantaheous 
response, in this case subsequent earthquakes should happen 
immediately after the event rather than being delayed for 
years. Aftershocks located at the end of the main shock 
fit this picture, however, the short time delay of the 
aftershock is most likely due to local inhomogeneities 
and creep type relaxation rather than large scale mantle 
relaxation. On the other hand, for the laterally inhomogeneous 
model G2, a significant portion of the perturbing stress in 
front of fault tip is accumulated years after the event. 

The perturbing stress levels off a few decades after the 

event. This accelerated stress accumulation in front of the 

region of a strike slip fault makes the chance of earthquake 
happening greater during the several decades after the event. 
Triggered events can then happen years after the triggering 
event, as long as the stress diffusion is reenforcing the 
prestress. The earthquake sequence after the 1939 North 
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Anatolian event (Toksfiz et al. , 1979) is consistent with this 
stress diffusion mechanism. 

It has been reported (Thatcher, 1975a) that the strain 
accumulation near the fault after the 1906 San Francisco 
event yas episodic: the average strain rate parallel to the 
fault was about 2.5 x 10“® per year for about 25 years after 

the event, and 0.6 x 10“® per year afterwards. For model Gl, 

! 

the average strain rate at a point 18 km away from the fault 

^ - . ■ . — g i 

xs about 0.1 X 10 per year, much less than the observed 
value. The relaxation does not have an episodic behavior 
either. This simple model does not explain the 1906 data. 

For the laterally inhomogeneous model G2 , the average strain 
rate at 18 km away from the“ fault is about 1.0 x 10 ° per 
year, similar to the observed values of the 1906 San Francisco 
earthquake. The tapering off of accelerated strain rate in 
about 30 years is also consistent with the San Francisco 
earthquake. The model is not a model of the San Francisco 
earthquake in detail, however, the results indicate that if 
the low viscosity extends to shallow depth under the fault, 
viscoelastic relaxation can contribute significantly to the 
time dependent deformation. 

Model G3 has the low viscosity channel to 40 km depth 
instead of 20 km in model G2. |The| results are intermediate 
between Gl and G2. Figure 13 compares the contours of 
vertical displacement change at 25 years after the event 
for model G2 and G3. The result for model G3 in the upper 
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right quadrant is a broad region of uplift, and a small 

' • -J 

region of slight subsidence near the fault tip. This 
indicates that time dependent vertical displacement is 
sensitive to the viscosity structure. Figure 14 compares 
the shear stress at selected positons vs. time for 
models G2 and G3. There is accelerated strain rate a few 
decades after the event for model G3, but the magnitude is 
smaller than that of model G2. 

i ' 

Moderate Size Strike Slip Earthquake Models 

Models and M 2 are models of moderate size earthquakes 

with buried faults. Model Ml is the fast relaxation model, 

IQ- 

with viscosity 10 poise in the low viscosity zone, whil 0 
model M 2 uses viscosity 10^® poise instead. The“ horizontal 

I 

displacements at selected locations on free surface are 

shown in Figures l5a and 15b. The results are generally a 

! 

gradual increase in magnitude. Figure 16 shows the 
instantaneous response of vertical deformation on free 
surface. In the upper right quadrant, there is a relatively 
broad region of uplift near the fault, and subsidence far 
away. This is expected for a fault with large fault depth : 
to length ratio (0.6 in this case). The vertical deformation 
changes at 9.5 years later are given in Figures 17a and 17b. , 

. -r ' " i 

For both models Ml and M2, in the upper right quadrant, ^ 

the vertical deformation is subsidence near the fault zone 
and uplift far away from the fault, similar to that of the 
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layered medium casei This iis probably because bhe width 
of the low viscosity channel is relatively large compared 
to the earthquake fault dimension. The magnitude of time 
dependent vertical deformation is small, generally less 
than 1 cm. ’ 

Figure 18 gives the time dependence of shear stress 

I • - . 

cT at selected locations at 7.5 km depth. We again see 
xy . 

the prestress is relieved along the fault, and reenforced 
in front of the fault.. The time dependence is quite 
different for the two models; significant changes are, 
completed within 5 years after the event for model Ml, 
while the relaxation is linear in modc‘1 M2 for the first 
decades. Though the magnitudes of quantities involved ib 
those moderate earthquakes are of marginal use with the 
precision of present day available geodetic data, the 
result may be useful in conjunction with more precise and 
fast geodetic measurements in the future to infer the 
detailed structure near a fault zone. 

IV. Discussion 

j 

In this study, we calculated time-dependent stress 
relaxation and deformation for large and moderate size 
earthquakes using different models . To compare these to 
observations we look at space-time migration of 
seismicity and to geodetic data. Migration of seismicity 
along plate boundaries has been observed in South America 
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(Kelleher, 1972), the West Pacific (Mogi, 1968), Alaska- 
Aleutian (Kelleher, 1970) , the San Andreas Fault zone 
(e.g. Wood and Allen, 1973; Lee et al., 1979), and the North 
Anatolian fault zone (e.g. Richter, 1958; Mogi, 1968; 

Savage, 1971; Dewey, 1976; Toksfiz et al., 1979). In the 
case of the North Anatolian faults there was a relatively 
quiescent period in seismicity prior to 1939. A bi- 
directional trend of seismic migration followed the great 
earthquake of 1939 (Toksflz et al., 1979). The occurrence of 
a sequence of events in rapid succession favors an accelerated 
stfess accumulation process. In all the model results in 
section III, there is an accelerated stress accumulation 
prqcess in the region in front of the fault tip. The | 

. ■ i i ; , ■ ; " ! ' ■ ' ' 

magnitude of this adcelerated stress accumulation is 
significant if a low viscosity zone under the fault extends 
to a depth of 40 km or less below the surface. 

The possible scenario for trie bccurrence of a sequence 
of earthquake follows. The initial earthquake happens when 


stress accvimulation exceeds the material strength. After 
this event, i trie stress accumulation accelerates in the 
region in front of the fault tip. Thfe next earthquake 
happens when the combinatidn Lof diffused stress and initial 


It , II 

stress exceeds the strength; in turn this triggered earthquake 
triggers the next events in an adjacent region. The process 
continues until the stress along the whole fault zone is 
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relieved. This episode is then followed by a slow stress 
accumulation stage and relative quiescence in seismicity. 
Aftershocks that immediately follow the earthquake are 
probably due to local stress adjustments. Creep resulting 
from stress changes at the immediate vicinity of the source 

I 

may result in relatively rapid stress adjustments in the 

source area. Migration of earthquakes in time may be 

■ - . I i 

related to viscoelastic relaxation and stress diffusion. 

Time-dependent horizontal and vertical motions after a 

! - - J _ ■ 'I ' ■ i 

strike slip extent strongly depend^^^^^^^^b^ ^ 
under and near the fault zone. Although horizontal^ 
displacements are much larger than vertical, in this study 
we found that the time dependent behavior of vertical 
displacement is very sensitive to lateral heterogeneities 
of viscosity distribution. The bulge or subsidence formed 
after a great earthquake is of meaeureable magnitude. Thus 
levelling, in addition to horizontal geodetic measurements, 
after a great strike-slip earthquake will reveal the structure 
near the fault zone. 

^ ■ 1 ■ 

Deep aseisnuc slip below the seismo-genic layer is 
suspected of playing an important role in earthquake mechanism 
Thatcher (1975a) reported accelerated strain rate decades 
prior to and after the 1906 San Francisco earthquake, and 
this can be explained by deep aseismic slip. Several studies 
have found that it is difficult to distinguish the effects 
of deep aseismic slip and viscoelastic relaxation from 
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geodetic measurements alone (Barker, 1976; Bundle and 
Jackson, 1977a, b; Savage and Prescott, 1978), However, 
this argument was based on calculations using laterally 
homogeneous layered models and comparisons of the surface 
deformation caused by viscoelastic relaxation and a certain 
stajtic dislocation at depth. As pointed out by Savage and 
Prescott (1978) , there are inherent difficulties in using 
a simple layered model to compare the results, because in 
such cases a distribution of slipj at depth can always 
duplicate the viscoelastic result'. This inhereht difficulty 

may not arise if the symmetry is broken by three dimensional 

^ ; i i . 

inhomogeneity. The relaxation is more intensive near a 
low viscosity zone, as we have shown in section III. 

Although present day available geodetic data are not very 
constraining, geodetic measurements in the future could 



A pre-seismic strain rate increase was reported in the 
case of the 1906 San Francisco earthquake (Thatcher, 1975a). 
This phenomenon cannot be caused by viscoelastic relaxation 
alorie, since relaxation is not a Spontaneous process. However 
if viscosity is indeed low near the ■ fault zoneV relaxation 
will add to the strain accumulation. A p^ossible w^ 
viscoelastic relaxation enters the accelerated strain rate 
process could be a renforcement effect culminati'ng in the 
fracture of the seismo-genic layer; the aseismic slip at 



depth increases the prestress on the locked fault, the effect 
of viscoelastic relaxation in general increases the 

i 

deformation, so the prestress near the fault is further 
increased. The magnitude and relative importance of this 
positive reenforcement again depends on the viscosity near 
the fault. 

In , conclusion, we have implemented a versatile scheme 
to model the time dependent behavior after earthquakes. 
Non-uniform fault slip and threie dimensional heterogeneities 
can be included in this scheme. The model results predict 
a stress diffusion phenomenon in front of fault tip after a 
strike slip event: if low viscosity extends to shallow depth 
near the fault zone, the shear stress in front of the fault 
tip will increase signif icant.ly with time. The time dependent 
deformation on free surface is more concentrated near the 
fault zone in that case than it is in the case of a laterally 
homogeneous layered structure. The time dependent behavior 
of vertical displacement near the fault may be completely 
altered by the presence of lateral inhomogeneities. 




APPENDIX A 

Time Dependent Finite Element Scheme Implementation 
^ We combined the time stepping approach of Zienkiewicz 
and Cormeau (1974) and frontal solution of Irons (1970) 
for our finite element calculations. The central process 
of the frontal solution approach is the familiar Gaussian 
elimination. However, a global displacement- force equation 
is eliminated i as soon as it has received all of the element 
and consistent nodal force contributions. The coefficients 
are moved to outside storage and other equations are updated 
accordingly. So only a small portion of the stiffness 
matrix will be in core at a time. This achieves the 
savings in core usage. Assembly and elimination processes 
are not separated. After completion of assembly and 
elimination, the coefficients are returned to core for , 
back substitution in the order exactly opposite to which 
they were saved. 

- i-r:-.. , ; j j . 

In the time stepping approach for time-dependent 
calculation the total strain e is divided into three parts; 

e = e^+ e'=P+ (A-1) 

where g® is the elastic strain, g° the initial strain arid 
the creep strain. Quite generally the constitutive 
law for creep strain rate and stress can be put into the form 

(A-2) 


- r a 
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where dot indicates time rate, r is a symmetric matrix (may 
depend on stress) and a is the stress. 

The virtual work principle is 

j 

(A-3) 

ISe^adQ, - /6u*^bdJi - /6u'^t d S = 0" ■ 

~ ~ ~ 

where b is the prescribed body force, t the prescribed boundary 

traction, u the displacement, 6 indicates variation and T indicates 
^ 

transpose. Integration is over the volume fi and traction 
boundary s^j respectively. 

i Let e=Lu, u = Na (A- 4) 

L is the operator relating strain and displacement, N the 
shape function and a the nodal displacements. Then equation 
vA3) becomes 

6a'^[/{L N)^a dfi - /n\ dfi - /N^ t dS] = 0 (A-5) 

^ ~ S'! "" 

or /b"^ a d^~ F = 0 

where B = L N, F = bdJ2+/N'^tdS 

^ Sg- 

Using a = D = D ( e - - e°) , equation (A5) can be put 

rw • -yj ~ ^ 

into standard form K a = V (A-6) 

where K = D B dJ^ - ; _ 

V = F + D e® d 0 + /b'^ D d 

~ ! Q~ " ^ ■' "J 

Equation (A6) is the set of linear equations to be solved 

using frontal solution, is obtained in time stepping 

fashion using equation = fa • it can be shown that this 

procedure can be applied to the general visco-plastic problem 

(Zienkiewicz and Cormeau, 1974) , including plasticity and 

creep problem as two extreme cases. The scheme can handle 
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npnlinear creep behavior as easy as linear material . However , 
nonlinear creep behavior depends on ambient 

stresses, since the strain rate depends on the sum of 
ambient and perturbing stresses. In this paper we have 
assumed that linearity holds and treated only the perturbation 
caused by the earthquakes. 
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Figure Captions 

Fig. 1. Fault for models Gl> G2 ^ and G3. 

i i 

(a) Schematic diagrams of fault, double hatched 
area has maximum fault slip, single hatched area 


has tapered fault slip. 

■' '^1 '5 

(b,) Fault slip along strike direption (Y direction) . 

(c) Fault slip vs. depth (Z direction). 

Fig. 2. Fault for models Ml and M2 

(a) Schematic diagram of fault, double hatched area has 
maximum fault slip, single hatched area has tapered 
fault slip. 

(b) Fault slip along strike direction (Y direction) . 

(c) Fault slip vs. depth (Z direction) . 

Fig. 3. Sectional views of viscosity distribution for models (a) G1 
(b) G2, (c) G3. Numbers with exponent are viscosity in poise. 

Fig. 4. Sectional view of viscosity distribution for model 

Ml. Model M2 has the same structure except that the 

• . 20 

vxscosity IS 10 poxse xn low vxscosxty zone. 

Fig. 5. Top view of the finite element grid. The three 

dimensional model is made of seven identical plane grids. 

Fig. 6. (a) Horizontal displacements on free surface due to 

strike slip event in model G1 at 0, 9.5, 25 and 49 
years after the event. The location of the fault is 
indicated by a thick line segment and sense of motion 
is indicated by a pair of arrows. 
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Dots are the locations where displacements are 
calculated. Displacement is indicated by a line . 
segment from the dot. A scale for the displacement 
(100 cm) and a scale for the map (400 km) are also 
shown in the figure. 

(b) The same plot for model G2. 

Fig. 7. (a) Horizontal displacement vs. time along the line 

perpendicular to the center of fault (x axis) on free 
surface for model Gl. The distance from the center of 
the fault for points 1, 2, 3, 4, 5, 6 and 7 are 35, 70, 
105, 140, 210 and 280 km. Schematic diagram of the 
locations are shown at the bottom of the graph. 

(b) The same plot for model G2. 

Fig. 8. (a) Contours of vertical displacement on free surface 

immediately after event for great earthquake models 
Gl, G2, and G3 . Broken lines (negative values) indicate 
subsidence; solid lines (positive values) indicate uplift. 
The location of fault is indicated by a thick line segment 
and a pair of arrows . The numbers near the contours are 
uplift or subsidence in mm. The tick marks on the 
frame are at half fault length interval (175 km). The 
elastic response is the same for model Gl, G2, and G3. 

(b) Contours of shear stress component Cj-y at 10 km 
depth inmiediately after the event for models Gl, G2, and 
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G3. The; nioitibers near the contours are stress in bars. 
Stress concentrations are placed near fault tips from 
intejrpblation, the actual grid resolution Is 1/10 of 
the fault length. ^ 

Fig. 9. (a) The vertical displacement on free surface at 

5 years after the event minus the elastic response for 
model Gl. Numbers are amount of uplift of subsidence 
in mm. Ticks are at half fault length intervals . 

(b) The same plot at 25 years after the event. 

Fig. 10. (a) Vertical displacement on free surface at 5 

years after the event minus the elastic response for 
model G2 . - 

! ' ! / 

(b) The same plot at 25 years after,.the bveiit. 

' , j ■ ’ 

Fig. 11. (a) Contours of shear stress component a^y at 10 km 

depth 25 years after the event for model Gl. 

(b) Same plot for model G2. r i = 

Fig. 12. (a) Shear stress component a^y vs. time j at selected 

locations for model Gl. The locatiohs are at 10 km depth. 
A schematic diagram indicating the horizontal positions 
relative to the fault is given to the right of the 
figure. The distances from the fault center along the . 
strike direction for points 1, 2, 3, 4> 5, 6,7, 8, ahd 
9 are 18,:, 53, 88, 123, 193, 228, 263, 307 and 385 km. 

The distance from the fault perpendicular to the strike 
direction is 35 km for point 9, 18 km for the rest of 
the points . 

(b) Same plot for model G2. 
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Pig. 13. (a) Contour of vertical displacement on free surface 

at 25 years after the event minus the elastic response 
for model G2 (Identical to figure 10b, included for 
comparison) . ^ , 

(b) The same plot for model G3. 

Fig. 14. 1 (a) Shear stress component cfj.y vs. time at 

selected locations for model G2. The locations are 

I • 

at 10 km depth. A schematic diagram indicating the 
horizontal positions relative to the fault is given to 
the right of the graph (identical to figure 12b, included 
for comparison) . 

(b) The same plot for model G3. 

Fig. 15. (a) Horizontal displacements on free surface due 

to strike slip event in model Ml at 0 , 1, 3 and 9.5 years 
after the event. 

(b) The same plot for model M2. 

1 - ■ 

Pig. 16. Vertical deformation on free surface immediately 

after the event for models Ml and M2. Ticks are at 
half fault length interval (10 km) . Numbers near the 
contours are uplift (positive values) or subsidence 
(negative values) in mm. 

Fig. 17. (a) The vertical displacement on free surface at 

9.5 years after the event minus the elastic response 
for model Ml. Numbers are uplift or subsidence in 
mm. Ticks are at half fault length interval (10 km) . 
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(b) The same plot for model M2. - ' 

Fig. 18. (a) Shear stress component vs. time at selected 

loca.tions for model Ml. The locations are at 7.5 km 

depth. A schematic diagram indicating the horizontal 

positions relative to the fault is given to the right ' 

of the figure* The distance from the fault center slbng ) 

i . : I 

the strike direction for points 1, 2, 3, 4, 5, 6, 7, 8, 9 

and 10 is 1.3, 3.8, 6.3, 11.3, 13.8, 16.3, 18.8, 21.9, 

27 . 5 , and 33.8 km. The distance from the fault perpendicular 

to the strike direction is 3.8 km for point 9 and 10, and 

1.3 km for the rest of the points. 

(b) The same plot for model M2. 
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PART 3: TECTONIC STRESS - MODELS AND MAGNITUDES 

Introduction 

Understanding the state of stress in the Earth's 
lithosphere is one of the paramount problems in Earth 
tectonics. The stress state is linked to causes - loading 
and unloading, heating and cooling, plate motions and 
driving forces ^ to consequences - creep deformation and 
seismic failure, and to rheology - the depth over which 
stress can be supported and the time dependence of material 
properties. None of the classes of links has been 
characterized in sufficiently quantitative detail to define 
the stress tensor in the lithosphere without ambiguity and 

I 

without a long inference chain involving poorly ' 

tested assumptions. This paper deals with one cause of 

stress in the lithosphere : the system of forces that 

! 

maintain plate motions. Specifically addressed are ways 
by which models of tectonic stress in the plates can be 
used to constrain the magnitude of regional deviatoric 
stress in the Earth's lithosphere. 

Global models of the intraplate deviatoric stress that 
arises from the driving and resistive forces controlling plate 
motion have been given by Solomon et al . [1975] and Richardson - 

et al. [1976, 1979]. |A principal objective of those studiesj 
has been to find those sets of forces that best match the body 
of iritraplate stress observations. The observations held to be 
most reliable for such comparisons are the indications 
of principal stress directions inferred from the mechanisms of 
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midplate earthquakes , from in situ stress measurements, and 
from the strikes of stress-sensitive geological features. 

While a comparison of model predictions and observations on 
the basis of principal stress orientations is straightforward 

I 

and serves as a useful test by which to reject possible force 
models, such an exercise does not directly address the absolute 
magnitude of intraplate deviatoric stresses, since all 
deviatoric stresses in a model can be multiplied by an 
arbitrary constant without changing the relative magnitudes 
or bhe orientations of the principal stresses. We show in this 
paper, however, that under certain conditions the body of 
data on intraplate stress orientations does constrain the 
magnitude of tectonic stresses. i 

It might be argued that stress magnitudes can in principle 
be measured by direct in situ techniques in sufficient 
locations to charapterize the stress field for length scales 
comparable to plate dimensions, thus obviating the need to 

4 

apply , indirect arguments to constrain tectonic stress 
magnitudes. This eventud.ity is doubtful for the near term, 

because of the difficulty of extrapolating near-surface measurements 
to greater depths in the lithosphere, and because furthej: advances 
in technology /will be necessary to conduct routinely measurements : . , 
of in situ stress over the large fraction of the Earthy s surface 
covered by"oceans. On the basis of available 'hydrofracture data, 
it may at least be concluded that deviatoric Stress magnitudes 

j . ' \ 

are on the order of several hundred bars to depths of several 
kilometers in a number of continental regions [ Haimson , 1977 ; 

McGarr and Gay, 1978]. 
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Thus the question remains ; given the large and growing 
body of data on the orientations of principal stresses within 
the plates, what information on the magnitudes of regional 
deviatoric stresses qan: be obtained from numerical models 

I . 

for tectonic plate stresses? We discuss in this paper two 
routes by which useful information on stress magnitudes can 
be derived; (1) For the driving force models that best fit 
the stress orientation data, if independent information on 
the magnitude of one or more of the forces in the system can 
be obtained, then the magnitudes of the total predicted 
stress field are constrained to comparable precision. The 
best fitting force models we have examined to date all involve 
a significant contribution from iridqe forces, the pushing 
forces that arise because of the elevated topography of ridge 
axes with respect to abyssal sea floor. Since ridges exert 
forces equivalent to compressive plate stresses of 200-300 
bars magnitude', this 'leads to the prediction that regional 

I 

deviatoric stresses are of this magnitude. (2) If in the: 
vicinity of a known local source of stress, the observations 
of stress orientations indicate comparable control by the local 
and regional stress field, then the magnitude of the regional 
field may be estimated. This line of argument holds special 
promise for oceanic intraplate regions where earthquakes 
have occurred in the vicinity of islands or large bathymetric 
features characterized by sufficiently good topographic and 
gravity data to model the associated local lithospheric stress. 


It should be mentioned that when direct in siiu 
measurements of stress magnitudes have high reliability, 
the I magnitude data can be used alongside the stress 
orientation data as a more powerful set of constraints on 
both regional and local forces on the lithosphere. 
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Stress Magnitudes and Global ! Plate Models 


The comparison of predicted and observed directions of 
principal intraplate stresses can be a sensitive test of possible 
sources of stress. As noted above, if such a comparison indicate? 
a significant contribution from a source of stress of known or : 
estimable magnitude, then a strong constraint on the general mag- 
nitude of deviatoric stress in the lithosphere on regional scales 
is obtained. In this section, we summarize our recent work 
[ Richardson et al . , 1979] on testing global models of intraplate; 
stress predicted by plate tectonic driving forces against observed 
directions of principal stresses, with particular emphasis on 

.... j - 

possible inferences on the magnitude of deviatoric stresses. 

Premises . That observations of principal stress directions 
in the plates can be used to constrain plate tectonic driving 

■ . I 

force models requires the adoption of three working premises; 

(i) that regionally consistent stress orientation fields exist 
for large fractions of the stable interiors of plates, (ii) that 


such stress fields are steady over time periods less than that 

g ^ : j 

(v 10 years) characterizing changes in plate motions; and (iii) 
that a recognizable portion of these stress fields is dominated 
by contributions from plate tectonic forces. 

The first premise has substantia] obserra tional support for most 
the plates [ Sykes and' .Sbar , 1974; Sbar and : Sykes , 1973; ^ Richardson ! 
et al . , 1979] ; see Figure 1. The second premise depends on the 
question of whether in plate interiors _the deformation and stress 
arising from past plate boundary slip superpose to produce steady 
motion and stress, or whether individual stress 'waves' from large 
earthquakes are discernible [e.g., Anderson , 1975] . This issue 
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may be resolved by ultra-precise geodetic measurements of short- 
term plate motions soon to be made [ Niell ^ al • / 1979; Smith 
et ^. / 1979; Bender et , 1979 . The third premise will be the 
most difficult to establish with certainty, but is a reasonable 
working hpyothesis for regions well removed from such other 
notable sources of stress as recent tectonic or thermal activity, 
recent topographic loading or unloading, and pronounced structural 
heterogeneities . 

Possible Driving Forces . We consider several simply parame- | 

teri zed driving and resistive forces as potential elements of a 
- ' i 1 

plate tectonic force model: plate bbundary forces at ridges, 

trenches, transform faults and zones of continent-continent i 

collision; and basal forces associated with viscous interaction ' "' 


between the lithosphere and the asthenosphere . While all of 
these forces contribute to lithospheric stress, it is important 
to recognize that potentially large stress contributions can also 
arise from lithospheric cooling [ Turcotte and Oxburgh , 1973] , 
latitudinal plate motion [ Turcotte and Oxburgh , 1973] , crustal 
thickness inhomogeneities [Artyushkov, 1973], lithospheric loading 
and unloading [ Walcott , 19 70^ Watts and Cochran , 1974; Haxby and 
Turcotte, 1976] and past tectonic events [ Swolf s et al . , 1974; 
Tullis , 1977] . In the interpretation of stress observations in terras 
of plate driving forces, care must be exercised to remove or to 
avoid where possible the effects of these additional 'sources of 
lithospheric stress. 

The compressive stress produced in oceanic plates by the 
elevation of mid-ocean ridges is the easiest to quantify among 
the set of possible driving and resistive forces, and is in the 
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range 200-300 bars [ Hales , 1969; Frank s 1972; McKenzie ^ 1972] . 

At subduction zones, the negative buoyancy of subducted litho- 
sphere is bapable of exerting an extensional force equivalent to several 
kilobars: stress on the adjacent plates [ McKenzie , 1969; Turcotte 
and Schubert , 1971] , but the greater fraction of available 
pulling force is counterbalanced by forces resisting descent of 

I 

the slab into the mantle [ Smith and Toksflz , 1972;! Forsyth and 

I 

Uyeda , 1975; Richter , 1977]. The net pull . by Slabs on the 
surface plates is uncertain but is considerably smaller than that 

' ' ■ i 

due to available negative buoyancy. At zones of continent- 
continent collision, the net force on the .adjacent plates' may be 
resistive (net compression) ,L=because of the contribution from the 
excess tdpography of the mountain belt marking the collision 
zone; the contribution from topography 'iri'ydlves' shear stresses 
of 200-300 bars for j the main boundary fault at the base Of the 
Himalayas [ Bird , 1978] . , The resistive torce at 
transform faults is uncertain [ Brune et al. , 

190; Brace and Byerlee , 1970J , but is not likely to be a major 
contributor to the plate driving mechanism on the basis ot the 
relat0ely small fraction of ; boundary taken up by transforms for 
most plates and the poor correlation of plate jspeeds with length 

of transform boundary [ Forsyth and Uyeda , 19 75; Aggarwal jr - 1978] . 

;■ . : . . . . r '■ 1. 

The viscous traction at the' base of the plates is less well 

!■ i : i i. ' ' . 

characterized than plate boundary forces and is uncertain both in 

I 

magnitude and in direction. The uncertainties are linkeid to 
questions of the radial scale for upper mantle convection, the 
planform for ' counter flow' to balance plate creation and destruc- 
tion, and the existence of a smaller secondary scale of astheno- 
spheric convection to transport heat [ Richter and Parsons , 1975; 
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McKenzie and Weiss , 1975; Harper ^ 1978; Chase , 1979 ; Hager ! 
and O'Connell , 1979] . Some simple forms for viscous drag 
are adopted as a basis for testing models, but the various 
potential complexities must be kept in mind. 

Stress Models . A variety of driving force models 
incorporating different relative amounts of boundary forces 
and basal tractions as described above have been tested against 
the observations of intraplate stress orientations. The 
lithosphere is modeled as a thin, spherical, elastic shell, 
and stresses are calculated from the imposed forces using 
the finite-element : analysis described by Richardson [1978] . 

The results of many models are given in Richardson et ail. 

j 

[1979], and only a summary of the results pertinent to the 
question of stress magnitudes will be given here. 

A summary of stress orientation data for intraplate 
regions is given in Figure 1. Most of the data come from 
the mechanisms of intraplate earthquakes; the re- 
mainder are from in situ measurements (see Richardson 
et al . [1979] for the original sources of the data shown). 

In compiling such a data-set, it is necessary to establish criteria 
for the selection qf-thoste data most appropriate for constraining the 
tectonic stress field. While such criteria are of 

I , 

necessity at least partly arbitraicy, our approach has been 
to ex'clude only those data very near {a.100 km distance or 
less) plate boundaries and those data likely reflecting 
unmodeled processes. Thus data from continental margins 
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have been excluded on the basis of possible contributions 
from sediment loading or thermal contraction - effects not 
modeled, and data are not used from regions of complex 
tectonics not likely to be a response solely to plate-scale 
forces (e.g., Alps,' Appalachians, and North America west of 
the Rockies) . 

Based on a comparison with the observed stress orientations 
in Figure 1, the predicted stresses are in best agreement with 
the observations when pushing forces at ridges are included in 
the driving force model and when the net pulling force due to 
subducted lithosphere is comparable in magnitude or is at 
most a few times larger than Other forces acting on. the plates. 

On the basis of intraplate stresses, therefore, resistive 
forces opposing the motion of the slab with respect to the 
mantle must nearly balance the negative buoyancy of the 
relatively cool, dense slab,! in agreement with similar ! 
conclusions derived from other considerations [Smith and 
Toksflz , 1972; F or sy th and Uyeda , 1975, Richter, 1977] . 

The maximum ratio of net slab pull to net ridge push is 

not sensitive to a dependence of net slab pull on subduction rate 

or to the inclusion of other forces in the system. Forces 
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resisting further convergence at coritinental collision zonejs 
along the Eurasian plate are important for intraplate stresses , 
and improve the fit to the data in Europe, Asia, and the Indian 
plate. Resistive viscous drag forces acting on the base of the 
plate in a direction opposite to "absolute" plate velpeity 
improve the fit to the intraplate stress field for several plates 
(e.g., Nazca, South America). The intraplate stress field is relatively 
insensitive to an increased drag coefficient beneath old 
oceanic lithosphere compared to young oceanic or continental 
lithosphere. Increasing by a factor of five or ten the drag 
coefficient beneath continental relative to oceanic lithosphere 
changes the calculated stresses only slightly and has little 

^ i ^ - r ;i ^ 

effect on the overall fit to observed stresses as long as some 
resistive drag acts beneath oceanic plates. 

Models in which drag forces drive (i.e., act parallel to 
"absolute" plate velocity) rather than resist plate motions are 
in poor agreement with the data. This poor agreement may depend 
on the oversimplified model of the adopted interaction between 
the plate and the asthenosphere. As noted above, the actual flow 
pattern in the mantle, including counter flow and possible multiple 
scales of convection, may be 'considerably more complicated than has 
been assumed in these models. 

Two models that provide reasonably gdod fits to a large 
fraction of the intraplate stress orientation data are shown in 
Figures 2 and 3. In Figure 2 are shown the predicted intraplate 
stresses for a model with the following forces; (i) a symmetric 
pushing force at ridges equivalent to a compressive stress of 100 
bars across a 100 km thick plate, (ii) a synuftetric pulling force 
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at trenches of the same absolute magnitude, (iii) a symmetric 
resistive force at continental collision zones of the same 
absolute 'magnitude and (iv) a drag stress ^ where v is absolute 
plate velocity in cm/yr and D is 0 . Ijbaf/cm/yr beneath oceans and 
0.6 bar/cm/yr beneath continents. Note that only the relative 
magnitudes of these forces are constrained by the stress orienta- 
tion data; their absolute magnitudes are proportional to an un- 
specified multiplicative constant. 

, The predicted directions of principal stresses for 'this force 

model are in good agreement with the data for eastern North 

1 : ■ ■ ■ : ■ i i 

America, Europe, Asia near the Himalayas,- and the Indian plate. i 

The fit to the data is good in South America, especially far from 
the trench, and in western Africa and is acceptable in most of the 
Pacific plate. The orientation ofi the calculated maximum compres- 
sive stress in the Nazca plate for the model is only in moderate 
agreement with the orientation inferred from the single fault 

I _ ; . ! . 

plane solution available. The fit to the data in the northern 

i ^ . 

Pacifid, eastern Asia, and east Africa: is rather poor. The fit 
to the data in the northern Pacific and eastern Asia could 
be improyed if subduction zone or drag forces werb decreased along 
the western Pacific plate margin or if slab forces were concentrated 
on the subducted plate. No attempt, however, has been made to vary 
plate boundary forces locally to match inferred stresses. If such 
an approacJi were adopted, most observed stresses could 
be matched but the solution for the driving mechanism would be 
unjustifiably arbitrary and non-unique. 

In Figure 3 are shown the intraplate stresses for a force 
model that takes the approach of Davies [1978 ] and Richardson 
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[1978 ] based on the assumption that drag balances the net torque 
on each plate due to boundary forces. The resulting drag thus 
varies from plate to plate and need not bear a simple relation- 
ship to relative plate motions, in contrast to the drag derived 
from absolute plate motion models consistent with known relative 
velocities [ Solomon and Sleep , 1974; Solomon et al. , 1975; 

i 

Minster et al. , 1974]. The force model includes: (i) a symmetric 

force at ridges equivalent to a compressive stress of 100 bars 
across a 100-km thick plate, (ii) a symmetric resistive force at conti- 
nental convergence 'zones of twice this magnitude, ,(i.ii) a pulling force 
at trenches, on the subducted platS: only, equivalent to an exten- 
sional stress of 100 bars across a 100-km thick plate, and (y) a 
viscous drag on each plate, due to the rotation of the plate with 

respect to the underlying mantLe_ ;(which may be moving), determined 

: ; i ■ ■ 

by balancing the total vector torque on the plate from boundary forces . 

The predicted stress directions for this model (Figure 3) 
agree very well with the data for several areas. In the North 
American and Nazca plates, the orientation of the maximum compres- 
sive stress is well matched by the model. The fit is almost as 

good in Europe and in Asia north of the Himalayas.; In the Indian 

' ! ^ ^ ^ 

plate, compressive stresses trend NW-SE in; continental India, in 

agreement with the data, but the fit is poorer in Australia." In 
South America, the maximum compressive stress trends E-W, I in only 
moderate agreement with the data^ In the Pacific and the; eastern 
part of the African plate the agreement with the data is poor. On 

the whole , this model provides a better f it i to continental than 

- ' i i ■ 

oceanic data. Comparison of Figures 2 and 3 suggests that any force 
pulling the overthrust plate toward the trench is probably lower in 
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magnitude than the net pull on the subducted plate. 

Discussion . From the standpoint of deviatoric stress 
magnitudes, the most important general conclusion from the 
modeling of plate tectonic stresses and the comparison with 

I 

intraplate stress orientation data is that ridge pushing forces 
are an imlport^n^ element of the set of driving forces for the 
models that provide the best fit to observations. The stresses 
that arise from ridge topography are 200-300 bars compression, 
as noted above. We are thus led to the conclusion that regional 
deviatoric stresses in plate interiors are of this same general j 
magnitude, or 200-300 barSi to within |a factor of perhaps 2 to 3. 

This conclusion should be tempered,: however; by several 
general observations on the results of the plate tectonic stress 
models. The models | represented in Figures 2 and 3/ though pro- 
viding good matches to the data for a number of regions with well 
characterized stresses, do not fit all of the data. Thus either 
there are simple models not tested that provide a better fit to 


the global data set than those: shown, or the stress observations 
are influenced by processes not included in the simpler models. 
Even if a model were obtained that fit all reliable observatiohs 


to within their estimated . errors, it is likely on the basis 
of models tested to date that this model would not be unique. 


Thus statements based on elements of best-fitting force models 
must be understood in recognition of this nonunigueness . 
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Stress Magnitudes and Local vs. Regional Stresses ' 

An alternative approach to constrain the magnitude of 

regional deviatoric stresses in the lithosphere from stress 
: ■ .1 
orientation data and plate tectonic models | is to find situations 

in which observed stress orientations are sensitive in approx- 
imately equal measure to a local stress field that may be readily 
quantified and to a regional stress field whose magnitude is to 
be determined. Such an approach holds high promise for constrain- 
ing the magnitudes of plate tectonic stresses in oceanic litho- 
sphere. 

Consider the effect of a volcanic load on oceanic lithosphere. 
Such a load leads to lithospheric flexure and to potentially large 
local: bending stresses. For a very large load, such as Hawaii, 
the local stresses may be in excess: of 1 kbar [ Walcott , 197:0;: 

Watts and Cochran , 1974] and may dominate the regional stress. I . 
That bending stresses may dominate regional stresses for Hawaii 
is supported by the report by Rogers and Endo : [1977] that greatest 
compressive stress axes from composite fault plane solutions for 
many mantle earthquakes beneath and near the island of Hawaii are 
radial with respect to the; island. 

I For loads appropriately smaller in magnitude than Hawaii, 
the local stresses should be comparable in magnitude to 
the regional stresses. Thus the mechanisms of earthquakes in the 
vicinity of such loads might be expected to indicate P and T axes 
which differ somewhat from regional trends but which are not pre- 
dictable simply from stress models for the local load only. For 
earthquakes near very small loads or distant from any pronounced 
topographic relief, the mechanisms should reflect the regional 
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stress field. 

As an illustration of this approach, consider the region 
near the Ninetyeast ridge in the central Indian Ocean. The j 
Ninetyeast ridge is a pronounced linear feature some 5000 km 
long and rising 1500-2000 m above the surrounding seafloor 
[e.gj .,5 Bowin , 1973]. The ridge is isostatically compensated 
except at short wavelengths [ Bowin , 1973; De trick and watts , 

1979] . Several large earthquakes have occurred in the Indian 
plate in the general vicinity during this century [ Sykes , 

1970; Stein and Qkal , lj978] . 

The orientation of principal stresses in the Indian plate 
may 'be estimated from the fault plane solutions Of intraplate 
earth(^uakes. Figure 4 shows the P axis orientations for all large 
earthquakes with: known focal mechanisms in the Indian plate near the 
Ninetyealst ridge. There is a strongly regionally consistent 
NW-SE trend to the . direction of inferred greatest compressive 
stress . 

Two aspects of this general consistency _are , h<itewor thy : 

(i) The P axes for strike-slip events on and: near the Ninetyeast 
ridge trend in general agreement with those f on thrust .events 
in the plate off the ridge. Thus while a zone of weakness 
associated With the ridge may control the type of faulting 
[ Stein and Qkal , 1978] , the inferred direction of maximum 
horizontal stress for Ninetyeast ridge events is still 
reliable. The data in Figure 4 are entirely consistent 
with a generally uniform stress field across the portion of 
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the Indian plate shown, with strike-slip rather than thrust 
motion the preferred fault type within weak zones in the 
lithosphere. (ii) The P axes for thrust events off the 
Ninetyeast ridge are not orthogonal to the strike of the 
ridge. Thus stresses associated with ridge topography do 
not dominate the local stress field. 

This second conclusion can be quantified to produce a 
constraint on the magnitude of the regional stress field. 
Adopting Bowin* s [1973] model for the isostatic compensation 
of the Ninetyeast ridge, the compressive force that the 
ridge exerts per unit length on the lithosphere beneath the 


adjacent abyssal plain may be estimated from equations (47- 
49) in Artyushkov [ 1973] : 

E = / (a - a ) dz 

XX zz 


( 1 ) 


where and a, 

z z 


(.14? + .067?^) X 10 bar-cm + Xj-idge 
are horizontal and vertical normal stress 
components ('\j principal stresses), ? is the height of the ridge 
(in km) with respect to the abyssal plain, is the value 

of E beneath the ridge, the integral is taken over the depth 
range of horizontal density variations, and the minus sign 
denotes a compressive force. ^ Notd^ that (1) includes the 
effects of topography and isostatic compensation only; the 
effects of viscous forces at the base of the plate and of 
thermal stress due to any differential cooling between the 
ridge and surrounding sea floor, for instance, are not included 
For ? = 1.5 to 2 km [ Bowin , 1973], (1) gi'^^ss E-E^j^^^^ = 


- (0.35 to 0.54) X 10^ bar-cm, or the equivalent of 70 
to 110 bars additional horizontal deviatoric stress over a 
50 km thick plate, for comparison, Artyushkov [1973] gives 

g 

“ 1.2 X 10 bar-cm and 240 bars compression for the force/ 
length and stress associated with spreading ridges. 

Thus the regional deviatoric stresses in the Indian plate 
(excluding the contribution from the Nihetyeast ridge) must 
be larger than ~100 bar in magnitude in order to account for 
the pattern of stress orientations ini Figure 4. This result 
provides only a lower bound on the magnitude of regional 

deviatoric stresses in one plate., but the result 

is at least consistent with the inference made above that 
regional stresses are similar in magnitude to the stresses 
produced by ridge forces, which are 3 + 1 times as large as 
the force exerted by Ninetyeast ridge topography. . 

A number of other oceanic intraplate earthquakes large 


enough so that their focal mechanisms are known have occurred 
in close proximity to prominent bathymetric features. 


Bergman arid Solomon [1980] have compiled a comprehensive 
catalog of 159 oceanic intraplate earthquakes and, for 
a representative s\abset of 83 epicentral regions, have 
assessed the degree of association of these earthquakes 
with large bathymetric features and with zones of expected 
lithospheric weakness (e.g., fracture zones). A histo- 
gram of the results (Figure 5) shows that the epicenters 
of at least 12 oceanic intraplate earthquakes with known 
focal mechanisms, and at least 11 additional earthquake 
localities with one or more events since 1964, 

are situated near features of pronounced seafloor topography. 
Several of these features involve lithospheric loads that 
should lead to bending stresses larger than the stresses 
indicated above for the Ninetyeast ridge. Thus it may be. 
possible by a combination of detailed stress models and 
careful source mechanisms to bound regional deviatoric 
stress magnitudes from both above and below usin^ this 
approach. 


Two potential drff iculties with this approach 'should, ^ 
however, be noted:; (i)- Many oceanidi intraplate earthquakes 
occur in or near such obvious zones, of weakness as fracture 
zones and volcanic areas. Over 70 percent of the oceanic 
earthquake epicentral regions in the listing of Bergman and 
Solomon [1980] are located in such areas (Figure 5) . Stress directions 


inferred from earthquake mechanisms for such events should be 
used only with caution in the absence of corroborative 
information friom events removed ftoin tHe weak zone (e.g.^ 

Figure 4). (ii) Bending stresses associated with lithospheric 
flexure, for given rheology and thermal structure, are 'extremely j 
sensitive to depth. Thus for an observation of ptress orientation from 
an earthquake mechanism to be a useful constraint on stress amplitude, 
the focal depth must be known with high precision, probably to 
within a few kilometers . 

Conclusions 

The global data on directions of principal stresses in 


plate interiors can serve as a test of possible plate tectonic 
force models. Such tests conducted to date favor force models 


in which ridge pushing forces play a Significant role. For 
such models, the general magnitude of regional deviator ic 
stresses is comparable to the 200-300 bars compressive stress 
exerted by spreading ; ridges . 


An alternative approach to estimating magnitudes of 
regional de^jiatoric stresses from stress ori|entations is to 
look for regions of local stress either demonstrably smaller 


than or larger than the regional stresses. The regional 
stresses in oceanic intraplate regions are larger than the 
'^'100 bar compression exerted by the Ninetyeast ridge and less 
•:han the bending stresses (>1 kbar) beneath Hawaii. 
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Figure 1. A summary of intraplate stress orientation data 
[ Richardson et al ., 1979], Filled circles denote fault 
plane solutions; arrows denote P and T axes^ where nearly 
horizontal. Filled circles without arrows denote thrust 
faults with poorly constrained P axes. Open circles 
represent in situ data; the line gives the direction of 
maximum horizontal compressive stress. 

Figure 2. Principal horizontal deviatoric stresses in the 

lithosphere for a model of plate driving forces (see text) . 
Principal stress axes without arrows and with arrows pointing 
outward denote deviatoric compression and tension, respec- 
tively. Relative magnitude of principal stresses is indicated 
by the length of stress axes. From Richardson et al . [1979]. 

Figure 3. Principal horizontal deviatoric stresses in the litho- 
sphere for an alternative driving force model in which basal 
shear balances the torque due to boundary forces for each 
plate (see text) . From Richardson et| al ,.,. [1979] . 

Figure 4. Summary of focal mechanisms for earthquakes in the Ninety- 
east ridge region of the Indian plate. Fault plane solutions 
are shown as equal area projections; compressional quadrants 
are shaded. Lines through each solution denote the orienta- 
tion of the P axis. Data are from Sykes [1970],. Fitch [1972], 
Sykes and Sbar [1974], Stein and Okal [1978], and Bergman and 


Solomon [1980] . 






Figure 5. Histograms of the number of oceanic intraplate earthquake 


epicentral regions in the catalog of Bergman and Solomon 


[1980] sorted by likelihood 


of association with either (left) 


pre-existing zones of weakness, such as fracture zones or 


volcanic seamount chains, or (right) with topographic relief 
that may provide a significant local source of stress. 
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